# Project: Don't Trust North Korea until It Sheds Its Skin
##         large North Korea's Responses to UN Human Rights Recommendations
# Author(s): Sanghoon Park & Hyunkyu Kim
# Time-log: 2024-Aug

# Load packages ---------------------

## Plot theme set
theme_clean <- function() {
  theme_minimal(base_family = "Barlow Semi Condensed") +
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(face = "bold", color = "black"),
          plot.subtitle = element_text(face = "bold", color = "black", size = rel(1.4)),
          axis.title.x = element_text(family = "Barlow Semi Condensed Medium", color = "black",
                                      margin = margin(t = 10, b = 10), size = rel(1.4)),
          axis.title.y = element_text(family = "Barlow Semi Condensed Medium", color = "black",
                                      margin = margin(r = 10, l = 10), size = rel(1.4)),
          axis.text = element_text(color = "black", size = rel(1.3)),
          strip.text = element_text(family = "Barlow Semi Condensed",
                                    face = "bold", size = rel(1), hjust = 0,
                                    color = "black"),
          strip.background = element_rect(fill = "white", color = NA),
          plot.caption = element_text(hjust = 0, color = "black"),
          legend.position = "bottom")
}
showtext::showtext_auto(T)
theme_set(ggpubr::theme_pubr())

# Import data -----------------------
## Sample data import --------
sample <- readRDS("Data/sample.RDS")

# DV -------------------------------------------------------------
sample |> 
  mutate(accept = if_else(response == "Supported", 1L, 0L),
         sensitivity = if_else(sensitivity > 1L, 1L, sensitivity)) -> 
  sample

sample |> dplyr::filter(year > 2009) |> 
  group_by(accept) |> count() |> drop_na() |> 
  ungroup() |> 
  mutate(total = sum(n),
         ratio = n/total,
         accept = factor(accept, 
                         levels = c(0, 1), label = c("Reject", "Accept"))) |> 
  ggplot(aes(x = accept, y = ratio, fill = accept, color = accept)) + 
  geom_col(show.legend = F) +
  labs(x = NULL, y = "Percentage\n") +
  geom_text(aes(label = paste0(round(ratio, 4)*100, "%", "\n(", n, ")")),
            position = position_stack(0.8), show.legend = F, 
            color = "white") +
  scale_color_manual(values = c(futurevisions::futurevisions("mars")[1],
                                futurevisions::futurevisions("mars")[3])) +
  scale_fill_manual(values = c(futurevisions::futurevisions("mars")[1],
                                futurevisions::futurevisions("mars")[3])) +
  scale_y_continuous(labels = scales::percent_format())

ggsave("Documents/Figures/Figure_DV.pdf", 
       width = 7, height = 4, dpi = 1600)

sample |> drop_na(REV_COW) -> sample

## Descriptive Statistics by DV -------

sample |> dplyr::filter(year > 2009) |> 
  dplyr::filter(accept == 1L) |> 
  dplyr::select(sensitivity, specificity, ideap_w_nk, ideap_w_ch, 
                REV_hrs, REV_dem, REV_gdppc, REV_gdpgrth, 
                l_aid_inflow) |> drop_na() |> 
  psych::describe()

sample |> dplyr::filter(year > 2009) |> 
  dplyr::filter(accept == 0L) |> 
  dplyr::select(sensitivity, specificity, ideap_w_nk, ideap_w_ch, 
                REV_hrs, REV_dem, REV_gdppc, REV_gdpgrth, 
                l_aid_inflow) |> drop_na() |> 
  psych::describe()

## Specificity and accept -----
table(sample$specificity[which(sample$year > 2009)])
sample |> dplyr::filter(year > 2009) |> 
  drop_na(accept, specificity) |> 
  group_by(specificity) |> 
  summarize(
    mean_accept = mean(accept, na.rm = T),
    sd_accept = sd(accept, na.rm = T),
    n = n(),
    se_accept = sd_accept/sqrt(n)
  ) -> specificity_accept

sample |> dplyr::filter(year > 2009) |> 
  drop_na(accept, specificity) |> 
  group_by(sensitivity) |> 
  summarize(
    mean_accept = mean(accept, na.rm = T),
    sd_accept = sd(accept, na.rm = T),
    n = n(),
    se_accept = sd_accept/sqrt(n)
  ) -> sensitivity_accept

sample |> mutate(extra = if_else(sensitivity == 0L & specificity == 0L,
                                 1L, 0L)) |> 
  drop_na(extra, accept) |> 
  group_by(extra) |> 
  summarize(
    mean_accept = mean(accept, na.rm = T),
    n = n()
  ) -> extra_accept

specificity_accept |> 
  ggplot(aes(x = specificity, y = mean_accept)) +
  geom_col(show.legend = F, fill = "#E7B800") +
  labs(x = NULL, subtitle = "Specificity") +
  geom_errorbar(aes(ymin = mean_accept - se_accept,
                    ymax = mean_accept + se_accept), width = 0.2) +
  geom_text(aes(label = paste0(round(mean_accept, 4)*100, "% (", 
                               n, ")")),
            position = position_stack(0.5),
            color = "black", size = rel(4)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(y = "Average Acceptance Rate") +
  scale_x_continuous(breaks = c(0, 1), labels = c("Non-specific", "Specific")) +
  theme(plot.subtitle = element_text(size = rel(1.5))) -> rec_specific_accept
  
sensitivity_accept |> 
  ggplot(aes(x = sensitivity, y = mean_accept)) +
  geom_col(show.legend = F, fill = "#FC4E07") +
  labs(x = NULL, subtitle = "Sensitivity") +
  geom_errorbar(aes(ymin = mean_accept - se_accept,
                    ymax = mean_accept + se_accept), width = 0.2) +
  geom_text(aes(label = paste0(round(mean_accept, 4)*100, "% (", 
                               n, ")")),
            position = position_stack(0.5),
            color = "black", size = rel(4)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(y = "Average Acceptance Rate") +
  scale_x_continuous(breaks = c(0, 1), labels = c("Non-sensitive", "Sensitive")) +
  theme(plot.subtitle = element_text(size = rel(1.5))) -> rec_sensitive_accept


rec_specific_accept + rec_sensitive_accept + 
  patchwork::plot_layout(ncol = 2, axes = "collect") 

ggsave("Documents/Figures/Figure_bivariate_accept_shaming.pdf", 
       width = 7, height = 5, dpi = 1600)
  
## Sensitivity and accept -----


sensitivity_accept |> 
  ggplot(aes(x = sensitivity, y = mean_accept)) +
  geom_col(show.legend = F, fill = "#FC4E07") +
  labs(x = NULL, subtitle = "Sensitivity") +
  geom_errorbar(aes(ymin = mean_accept - se_accept,
                    ymax = mean_accept + se_accept), width = 0.2) +
  geom_text(aes(label = paste0(round(mean_accept, 4)*100, "% (", 
                               n, ")")),
            position = position_stack(0.5),
            color = "black", size = rel(4)) +
  scale_y_continuous(labels = scales::percent_format()) +
  labs(y = "Average Acceptance Rate") +
  scale_x_continuous(breaks = c(0, 1), labels = c("Non-sensitive", "Sensitive")) +
  theme(plot.subtitle = element_text(size = rel(1.5)))

# Predictors -------

sample |> drop_na(REV_COW) |> 
  dplyr::filter(year > 2009) |> 
  dplyr::select(sensitivity, specificity) |> 
  tidyr::pivot_longer(
    cols = c("sensitivity", "specificity"),
    values_to = "score"
  ) |> group_by(name, score) |> 
  mutate(n = n()) |> drop_na() |> 
  group_by(name) |> 
  summarize(mean = mean(score, na.rm = T),
            n = mean(n)) |> 
  ggplot(aes(x = name, y = mean)) + geom_col()
library(janitor)
sample |> drop_na(REV_COW) |> 
  dplyr::filter(year > 2009) |> 
  mutate(#sensitivity = factor(sensitivity, 
        #                      levels = c(0, 1), labels = c("Non-sensitive", "Sensitive")),
         #specificity = factor(specificity,
        #                      levels = c(0, 1), labels = c("Non-specific", "Specific")),
         nothing = case_when(sensitivity == 1L & 
                               specificity == 0L ~ "Only Sensitive",
                             sensitivity == 0L & 
                               specificity == 1L ~ "Only Specific",
                             sensitivity == 1L & 
                               specificity == 1L ~ "Sensitive and Specific",
                             T ~ "Others")) |> 
  drop_na(nothing) |> 
  janitor::tabyl(nothing, show_missing_levels = F) |> 
  adorn_totals(c("row")) |> 
  adorn_pct_formatting(digits = 2) -> tbl1
tbl1 |> 
  kableExtra::kbl("latex", booktabs = T) |> 
  kableExtra::kable_paper(full_width = F)


library(rstatix);library(ggpubr)
sample |> dplyr::filter(year > 2009) |> 
  rstatix::t_test(ideap_w_nk ~ accept, ref.group = 1) ->
  t_test_accept_ideap_nk

t_test_accept_ideap_nk

my_comparisons <- list( c("1", "0"))

sample |> dplyr::filter(year > 2009) |> drop_na() |> 
  ggboxplot(x = "accept", y = "ideap_w_nk",
            color = "accept", palette =c(#"#00AFBB",
              "#E7B800", "#FC4E07"),
            add = "jitter", shape = "accept", legend = "none") + 
  scale_x_discrete(
    labels = c("Reject", "Accept")) +
  scale_y_continuous(limits = c(0, 6)) +
  labs(x = NULL, y = "Average Values\n",
       subtitle = "Geopolitical Affinity") +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  stat_compare_means(label.y = 5.5, label.x = 1.2) -> t_test_p1                  

ggsave("Documents/Figures/Figure_bivariate_accept_ideap_nk.pdf", 
       t_test_p1,
       width = 3.5, height = 4.5, dpi = 1600)

sample |> dplyr::filter(year > 2009) |> 
  rstatix::t_test(ideap_w_ch ~ accept, ref.group = 1) ->
  t_test_accept_ideap_ch

my_comparisons <- list( c("1", "0"))

sample |> dplyr::filter(year > 2009) |> drop_na() |> 
  ggboxplot(x = "accept", y = "ideap_w_ch",
            color = "accept", palette =c(#"#00AFBB",
              "#E7B800", "#FC4E07"),
            add = "jitter", shape = "accept", legend = "none") + 
  scale_x_discrete(
    labels = c("Reject", "Accept")) +
  scale_y_continuous(limits = c(0, 6)) +
  labs(x = NULL, y = "Geopolitical Affinity\n",
       subtitle = "With China") +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  stat_compare_means(label.y = 4.2, label.x = 1.2) -> t_test_p2                  

ggsave("Documents/Figures/Figure_bivariate_accept_ideap_ch.pdf", 
       t_test_p2,
       width = 3.5, height = 4.5, dpi = 1600)

t_test_p1 + t_test_p2 + 
  patchwork::plot_layout(ncol = 2, axes = "collect") 

ggsave("Documents/Figures/Figure_bivariate_accept_ideap.pdf", 
       width = 7, height = 4.5, dpi = 1600)

sample |> mutate(
  shared_gdp = REV_gdp/(SUR_gdp + REV_gdp)
) -> sample
  
sample |> dplyr::filter(year > 2009) |> 
  rstatix::t_test(shared_gdp ~ accept, ref.group = 1) ->
  t_test_accept_shared_gdp

t_test_accept_shared_gdp

sample |> dplyr::filter(year > 2009) |> 
  rstatix::t_test(REV_share_gdp_ch ~ accept, ref.group = 1) ->
  t_test_accept_shared_gdp_ch

my_comparisons <- list( c("1", "0"))

sample |> dplyr::filter(year > 2009) |> 
  drop_na(accept, shared_gdp) |> 
  ggboxplot(x = "accept", y = "shared_gdp",
            color = "accept", palette =c(#"#00AFBB",
              "#E7B800", "#FC4E07"),
            add = "jitter", shape = "accept", 
            legend = "none") + 
  scale_x_discrete(
    labels = c("Reject", "Accept")) +
  labs(x = NULL, y = "Average Values\n",
       subtitle = "Reviewer's GDP / Shared GDP") +
  scale_y_continuous(labels = scales::percent_format(),
                     breaks = c(seq(0, 1.2, 0.2)),
                     limits = c(0, 1.2)) +
  stat_compare_means(comparisons = my_comparisons,
                     label = "p.signif") +
  stat_compare_means(label.y = 1.1, label.x = 1.2) -> t_test_shared_nk

t_test_shared_nk

sample |> dplyr::filter(year > 2009) |> 
  drop_na(accept, REV_share_gdp_ch) |> 
  ggboxplot(x = "accept", y = "REV_share_gdp_ch",
            color = "accept", palette =c(#"#00AFBB",
              "#E7B800", "#FC4E07"),
            add = "jitter", shape = "accept", 
            legend = "none") + 
  scale_x_discrete(
    labels = c("Reject", "Accept")) +
  labs(x = NULL, y = "Reviewer's GDP / Shared GDP\n",
       subtitle = "With China") +
  scale_y_continuous(labels = scales::percent_format(),
                     breaks = c(seq(0, 1.2, 0.2)),
                     limits = c(0, 1.2)) +
  stat_compare_means(comparisons = my_comparisons,
                     label = "p.signif") +
  stat_compare_means(label.y = 0.65, label.x = 1.2) -> t_test_shared_ch

t_test_shared_ch

t_test_shared_nk + t_test_shared_ch + 
  patchwork::plot_layout(ncol = 2, axes = "collect") 

ggsave("Documents/Figures/Figure_bivariate_accept_economictie.pdf", 
       width = 7, height = 4.5, dpi = 1600)

library(patchwork)
t_test_p1 + t_test_shared_nk + 
  patchwork::plot_layout(ncol = 2, axes = "collect") 
ggsave("Documents/Figures/Figure_bivariate_affinity.pdf", 
       width = 7, height = 4.5, dpi = 1600)

t_test_p1
ggsave("Documents/Figures/Figure_bivariate_affinity.pdf", 
       width = 6, height = 4.5, dpi = 1600)


sample |> dplyr::filter(year > 2009) |> drop_na() |> 
  group_by(accept) |> summarize(
    mean = mean(REV_share_gdp_ch, na.rm = T)
  )

sample |> dplyr::filter(year > 2009) |> 
  rstatix::t_test(l_aid_inflow ~ accept, ref.group = 1) ->
  t_test_accept_aidfromch

t_test_accept_aidfromch

my_comparisons <- list( c("1", "0"))

sample |> dplyr::filter(year > 2009) |> 
  drop_na(accept, lnlaidinflow) |> 
  ggboxplot(x = "accept", y = "lnlaidinflow",
            color = "accept", palette =c(#"#00AFBB",
              "#E7B800", "#FC4E07"),
            add = "jitter", shape = "accept", legend = "none") + 
  scale_x_discrete(
    labels = c("Reject", "Accept")) +
  scale_y_continuous(limits = c(18.2, 19.2)) +
  labs(x = NULL, y = "Logged Aid Inflow from China\n") +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  stat_compare_means(label.y = 19.5, label.x = 1.2) -> t_test_aidch

t_test_aidch

sample |> dplyr::filter(year > 2009) |> 
  drop_na(accept, lnlaidinflow) |> 
  mutate(accept = if_else(accept == 1L, "Accepted", "Rejected")) |> 
  ggplot(aes(lnlaidinflow, color = accept, 
             fill = accept, group = accept)) + 
  geom_density(alpha = 0.3, show.legend = F) +
  #geom_histogram(aes(y = ..density..)) +
  scale_fill_manual(values = c("#E7B800", "#FC4E07")) +
  scale_color_manual(values = c("#E7B800", "#FC4E07")) +
  facet_wrap(~accept, ncol = 1)
  
  ggdensity(x = "lnlaidinflow",
            group = "accept",
            add = "mean", rug = TRUE,
            #color = "accept", fill = "accept",
            palette = c("#00AFBB", "#E7B800"))

ggsave("Documents/Figures/Figure_bivariate_accept_ideap_nk.pdf", 
       t_test_p1,
       width = 3.5, height = 4.5, dpi = 1600)

sample |> dplyr::filter(year > 2009) |> 
  drop_na(accept, REV_hrs) |> 
  ggboxplot(x = "accept", y = "REV_hrs",
            color = "accept", palette =c(#"#00AFBB",
              "#E7B800", "#FC4E07"),
            add = "jitter", shape = "accept", legend = "none") + 
  scale_x_discrete(
    labels = c("Reject", "Accept")) +
  scale_y_continuous(limits = c(0, 6)) +
  labs(x = NULL, y = "Human Rights Scores\n") +
  stat_compare_means(comparisons = my_comparisons, label = "p.signif") +
  stat_compare_means(label.y = 5.5, label.x = 1.2) -> t_test_hrs

t_test_hrs

ggsave("Documents/Figures/Figure_bivariate_accept_hrs.pdf", 
       t_test_hrs,
       width = 6, height = 4.5, dpi = 1600)


# Models ---------------------------------------------------------
## Benchmark models ----------------------------------------------

glm(accept ~ sensitivity, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench1

glm(accept ~ specificity, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench2

glm(accept ~ ideap_w_nk, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench3

glm(accept ~ ideap_w_ch, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench4

glm(accept ~ REV_hrs, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench5

glm(accept ~ shared_gdp, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench6

glm(accept ~ REV_share_gdp_ch, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench7

glm(accept ~ REV_gdpgrth, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench8

glm(accept ~ lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> bench9

texreg::texreg(list(bench1, bench2, bench3,
                    bench5, bench6, bench8, bench9),
                  custom.coef.names = c(
                    "(Intercept)"  = "(Constants)",
                    "sesitivity"   = "Sensitivity",
                    "specificity"  = "Specificity",
                    "ideap_w_nk"   = "Geopolitical Affinity with DPRK",
                    #"ideap_w_ch"   = "Geopolitical Affinity with China",
                    "REV_hrs"      = "Reviewer's Human Rights Score",
                    "shared_gdp"  = "Shared GDP of DPRK with Reviewers",
                    #"REV_share_gdp_ch"  = "Shared GDP of China with Reviewers",
                    "REV_gdpgrth"  = "Reviewer's Annual GDP Growth",
                    "lnlaidinflow" = "Logged China Aid to Reviewer"),
                  reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 1))

## Simple models ---------------------------------------------------

glm(accept ~ sensitivity + specificity, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple1

glm(accept ~ sensitivity + specificity + REV_hrs, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple2

glm(accept ~ sensitivity + specificity + ideap_w_nk, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple3

glm(accept ~ sensitivity + specificity + 
      ideap_w_nk + REV_hrs, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple4

glm(accept ~ sensitivity + specificity + 
      ideap_w_ch + REV_hrs, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple5

glm(accept ~ sensitivity + specificity + shared_gdp, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple6

glm(accept ~ sensitivity + specificity + REV_share_gdp_ch, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple7

glm(accept ~ sensitivity + specificity + REV_gdpgrth, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple8

glm(accept ~ sensitivity + specificity + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple9

glm(accept ~ sensitivity + specificity + 
      shared_gdp + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple10

glm(accept ~ sensitivity + specificity + 
      REV_share_gdp_ch + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple11

glm(accept ~ sensitivity + specificity + 
      ideap_w_nk + REV_hrs +
      shared_gdp + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple12

glm(accept ~ sensitivity + specificity + 
      ideap_w_ch + REV_hrs +
      REV_share_gdp_ch + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> simple13

texreg::texreg(list(simple1, simple2, simple3, simple4, simple6,
                    simple8, simple9, simple10, simple12),
               custom.coef.names = c(
                 "(Intercept)" = "(Constants)",
                 "sesitivity" = "Sensitivity",
                 "specificity" = "Specificity",
                 "REV_hrs" = "Reviewer's Human Rights Score",
                 "ideap_w_nk" = "Geopolitical Affinity with DPRK",
                 #"ideap_w_ch" = "Geopolitical Affinity with China",
                 "shared_gdp"  = "Shared GDP of DPRK with Reviewers",
                 #"REV_share_gdp_ch"  = "Shared GDP of China with Reviewers",
                 "REV_gdpgrth" = "Reviewer's Annual GDP Growth",
                 "lnlaidinflow" = "Logged China Aid to Reviewer"),
               reorder.coef = c(2, 3, 4, 5, 6, 7, 8, 1))

## Conditional models ----------------------------------------------

glm(accept ~ ideap_w_nk*(sensitivity + specificity) + REV_hrs, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> inter1

glm(accept ~ ideap_w_ch*(sensitivity + specificity) + REV_hrs, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> inter2

glm(accept ~ shared_gdp*(sensitivity + specificity) + ideap_w_nk, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> inter3

glm(accept ~ REV_share_gdp_ch*(sensitivity + specificity) + ideap_w_ch, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> inter4

glm(accept ~ REV_gdpgrth*(sensitivity + specificity) + ideap_w_nk, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> inter5

glm(accept ~ REV_gdpgrth*(sensitivity + specificity) + ideap_w_ch, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> inter6

texreg::texreg(list(inter1, inter2, inter3, inter4, inter5, inter6),
               custom.coef.names = c(
                 "(Intercept)" = "(Constants)",
                 "ideap_w_nk" = "Geopolitical Affinity with DPRK",
                 "sesitivity" = "Sensitivity",
                 "specificity" = "Specificity",
                 "REV_hrs" = "Reviewer's Human Rights Score",
                 "ideap_w_nk:sensitivity" = "Geopolitical Affinity with DPRK$\\times$Sensitivity",
                 "ideap_w_nk:specificity" = "Geopolitical Affinity with DPRK$\\times$Specificity",
                 "ideap_w_ch" = "Geopolitical Affinity with China",
                 "ideap_w_ch:sensitivity" = "Geopolitical Affinity with China$\\times$Sensitivity",
                 "ideap_w_ch:specificity" = "Geopolitical Affinity with China$\\times$Specificity",
                 "shared_gdp"  = "Shared GDP of DPRK with Reviewers",
                 "shared_gdp:sensitivity" = "Shared GDP of DPRK with Reviewers$\\times$Sensitivity",
                 "shared_gdp:specificity" = "Shared GDP of DPRK with Reviewers$\\times$Specificity",
                 "REV_share_gdp_ch"  = "Shared GDP of China with Reviewers",
                 "REV_share_gdp_ch:sensitivity" = "Shared GDP of China with Reviewers$\\times$Sensitivity",
                 "REV_share_gdp_ch:specificity" = "Shared GDP of China with Reviewers$\\times$Specificity",
                 "REV_gdpgrth" = "Reviewer's Annual GDP Growth",
                 "REV_gdpgrth:sensitivity" = "Reviewer's Annual GDP Growth$\\times$Sensitivity",
                 "REV_gdpgrth:specificity" = "Reviewer's Annual GDP Growth$\\times$Specificity"),
               reorder.coef = c(3, 4, 2, 6, 7, 8, 9, 10, 5, 11, 12, 13, 14, 15, 16, 17, 18, 19, 1))

## Full models only with DPRK ------------------

glm(accept ~ sensitivity + specificity + 
      ideap_w_nk + REV_hrs + REV_boix + 
      shared_gdp + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> full_nk1

glm(accept ~ ideap_w_nk*sensitivity + specificity + 
      ideap_w_nk + REV_hrs + REV_boix + 
      shared_gdp + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> full_nk2

glm(accept ~ sensitivity + ideap_w_nk*specificity + 
      ideap_w_nk + REV_hrs + REV_boix + 
      shared_gdp + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> full_nk3

glm(accept ~ ideap_w_nk*(sensitivity + specificity) + 
      ideap_w_nk + REV_hrs + REV_boix + 
      shared_gdp + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> full_nk4

sample |> mutate(sensitivity = factor(sensitivity, levels = c(0, 1), labels = c("Sensitive", "Non-sensitive")),
                 specificity = factor(specificity, levels = c(0, 1), labels = c("Specific", "Non-specific"))) ->
  sample

ggeffects::ggpredict(full_nk2, terms = c("ideap_w_nk", "sensitivity")) |> plot()
ggeffects::ggpredict(full_nk3, terms = c("ideap_w_nk", "specificity")) |> plot()

marginaleffects::plot_slopes(full_nk2, variables = "sensitivity", condition = "ideap_w_nk")
marginaleffects::plot_slopes(full_nk3, variables = "specificity", condition = "ideap_w_nk")


texreg::texreg(list(
  full_nk1, full_nk2, full_nk3, full_nk4),
  custom.coef.names = c(
    "(Intercept)" = "(Constants)",
    "sesitivity" = "Sensitivity",
    "specificity" = "Specificity",
    "ideap_w_nk" = "Geopolitical Affinity with DPRK",
    "REV_hrs" = "Reviewer's Human Rights Score",
    "REV_boix" = "Reviewer: Democracy",
    "shared_gdp"  = "Shared GDP of DPRK with Reviewers",
    "REV_gdpgrth" = "Reviewer's Annual GDP Growth",
    "lnlaidinflow" = "Logged China Aid to Reviewer",
    "ideap_w_nk:sensitivity" = "Geopolitical Affinity with DPRK$\\times$Sensitivity",
    "ideap_w_nk:specificity" = "Geopolitical Affinity with DPRK$\\times$Specificity"),
  reorder.coef = c(2, 3, 4, 10, 11, 5, 6, 7, 8, 9, 1))

## Full models only with CHN ------------------

glm(accept ~ sensitivity + specificity + 
      ideap_w_ch + REV_hrs + REV_boix + 
      REV_share_gdp_ch + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> full_ch1

glm(accept ~ ideap_w_ch*sensitivity + specificity + 
      ideap_w_ch + REV_hrs + REV_boix + 
      REV_share_gdp_ch + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> full_ch2

glm(accept ~ sensitivity + ideap_w_ch*specificity + 
      ideap_w_ch + REV_hrs + REV_boix + 
      REV_share_gdp_ch + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> full_ch3

glm(accept ~ ideap_w_ch*(sensitivity + specificity) + 
      ideap_w_ch + REV_hrs + REV_boix + 
      REV_share_gdp_ch + REV_gdpgrth + lnlaidinflow, 
    family = "binomial",
    data = sample |> dplyr::filter(year > 2009)) -> full_ch4


texreg::texreg(list(
  full_ch1, full_ch2, full_ch3, full_ch4),
custom.coef.names = c(
  "(Intercept)" = "(Constants)",
  "sesitivity" = "Sensitivity",
  "specificity" = "Specificity",
  "ideap_w_ch" = "Geopolitical Affinity with China",
  "REV_hrs" = "Reviewer's Human Rights Score",
  "REV_boix" = "Reviewer: Democracy",
  "REV_share_gdp_ch"  = "Shared GDP of China with Reviewers",
  "REV_gdpgrth" = "Reviewer's Annual GDP Growth",
  "lnlaidinflow" = "Logged China Aid to Reviewer",
  "ideap_w_ch:sensitivity" = "Geopolitical Affinity with China$\\times$Sensitivity",
  "ideap_w_ch:specificity" = "Geopolitical Affinity with China$\\times$Specificity"),
reorder.coef = c(2, 3, 4, 10, 11, 5, 6, 7, 8, 9, 1))


